AEI-2007-119 



Best network chirplet-chain: Near-optimal coherent detection of unmodeled 
gravitation wave chirps with a network of detectors 

Archana PaQ 

Max-Planck Institut fur Gravitationsphysik, Am Muhlenberg 1, 14476 Potsdam, Germany 

Eric Chassande-Mottirfj] 

CNRS, AstroParticule et Cosmologie, 10, rue Alice Domon et Leonie Duquet, 75205 Paris Cedex 13, France and 
Observatoire de la Cote d'Azur, Bd de I'Observatoire, BP 4229, 06304 Nice, France 

Olivier Rabast^D 

CNRS, AstroParticule et Cosmologie, 10, rue Alice Domon et Leonie Duquet, 15205 Paris Cedex 13, France 

(Dated: February 5, 2008) 

The searches of impulsive gravitational waves (GW) in the data of the ground-based interfer- 
ometers focus essentially on two types of waveforms: short unmodeled bursts from supernova core 
collapses and frequency modulated signals (or chirps) from inspiralling compact binaries. There is 
room for other types of searches based on different models. Our objective is to fill this gap. More 
specifically, we are interested in GW chirps "in general", i.e., with an arbitrary phase/frequency vs. 
time evolution. These unmodeled GW chirps may be considered as the generic signature of orbiting 
or spinning sources. We expect the quasi-periodic nature of the waveform to be preserved indepen- 
dently of the physics which governs the source motion. Several methods have been introduced to 
address the detection of unmodeled chirps using the data of a single detector. Those include the 
best chirplet chain (BCC) algorithm introduced by the authors. In the next years, several detec- 
tors will be in operation. Improvements can be expected from the joint observation of a GW by 
multiple detectors and the coherent analysis of their data namely, a larger sight horizon and the 
more accurate estimation of the source location and the wave polarization angles. Here, we present 
an extension of the BCC search to the multiple detector case. This work is based on the coherent 
analysis scheme proposed in the detection of inspiralling binary chirps. We revisit the derivation 
of the optimal statistic with a new formalism which allows the adaptation to the detection of un- 
modeled chirps. The method amounts to searching for salient paths in the combined time-frequency 
representation of two synthetic streams. The latter are time-series which combine the data from 
each detector linearly in such a way that all the GW signatures received are added constructively. 
We give a proof of principle for the full sky blind search in a simplified situation which shows that 
the joint estimation of the source sky location and chirp frequency is possible. 

PACS numbers: 04.80.Nn, 07.05.Kf, 95.55.Ym 



I. SUMMARY 

A large effort is underway to analyze the scientific 
data acquired jointly by the long-baseline interfcromet- 
ric gravitational wave (GW) detectors GEO 600, LIGO, 
TAMA and Virgo pQ. In this paper, we contribute to the 
methodologies employed for this analysis, and in partic- 
ular for the detection of impulsive GW signals. 

The current GW data analysis effort is targeted on two 
types of impulsive GWs. A first target is poorly known 
short bursts of GWs with a duration in the hundredth of 
a millisecond range. The astrophysically known sources 
of such GW bursts are supernovae core collapses (or other 
similar cataclysmic phenomenon). The second target is 
frequency modulated signals or chirps radiated by inspi- 
ralling binaries of compact objects (either neutron stars 
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(NS) or black holes (BH)). These chirp waveforms are 
well modeled and expected to last for few seconds to few 
minutes in the detector bandwidth. Our objective is to 
enlarge the signal range of impulsive GWs under consid- 
eration and to fill this gap between these two types. More 
specifically, we are interested in the detection of unmod- 
eled GW chirps which last from few tens of milliseconds 
to few seconds in the detector bandwidth. We shall de- 
tail in the next section the astrophysical motivation for 
considering this kind of GWs. 

Simultaneous observation and analysis of the jointly 
observed data from different GW detectors has definite 
obvious benefits. First and foremost, a GW detection 
can get confirmed or vetoed out with such a joint ob- 
servation. Further, the detector response depends on the 
position and orientation of the source and polarization of 
the wave. For this reason, the joint observation by mul- 
tiple detectors gives access to physical parameters such 
as source location and polarization. The use of multiple 
detectors also allows to enlarge the observational horizon 
and sky coverage. 

Till date, several methods have been proposed and im- 
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plemented to detect unmodeled GW chirps using a sin- 
gle detector which include the signal track search (STS) 
[2], the chirplet track search [3] and the best chirplct 
chain (BCC) search proposed by the authors [I]. How- 
ever, none of the above addresses the multiple detector 
case. This requires the designing of specific algorithms 
which are able to combine the information received by 
the different detectors. 

In practice, there are two approaches adopted to carry 
out network analysis from many detectors, coincidence 
and coherent approach. In coincidence approach, the 
data from each detector is processed independently and 
only coincident trigger events (in the arrival time and 
the parameter values) are retained and compared. On 
the other hand, in the coherent approach, the network 
as a whole is treated as a single "sensor" : the data from 
various detectors is analyzed jointly and combined into 
a single network statistic which is tested for detection. 
In the literature, it has been shown that the coherent 
approach performs better than the coincidence approach 
for GW short bursts [5] as well as GW chirps from coa- 
lescing binaries [B]. Indeed, the signal phase information 
is preserved with the coherent approach, whereas it is 
not with the coincidence approach which therefore leads 
to some information loss. Thus, we wish to follow the 
coherent approach for our analysis. 

Another reason for this choice is that the coincidence 
method is not adequate for unmodeled chirps. A large 
number of parameters (in the same order as the num- 
ber of signal samples) is needed to characterize their fre- 
quency evolution. A coincident detection occurs when 
the parameter estimates obtained from the analysis of 
the individual detector data match. Because of the noise 
fluctuation, the occurrence of such a coincidence is very 
unlikely when the number of parameters is large, unless 
the incoming GW has very large amplitude. In this arti- 
cle, we therefore adopt the coherent method and propose 
the coherent extension of the BCC algorithm. 

Coherent schemes have been already developed for the 
detection of inspiralling binary chirps [7J [5]. Here, we 
revisit the work presented in [7J with a new formalism. 
Comments in footnotes link the results presented here 
with the ones of [7J. We show that the new formalism 
presented here helps to understand the geometry of the 
problem and it is simple to establish connections with 
earlier works. 

The outline of the paper goes as follows. In Sec. |H] 
we present and motivate our model of an arbitrary GW 
chirp. In Sec. |III[ we describe the response of the detec- 
tor network to an incoming GW chirp. Further, we show 
that the linear component of the signal model (parame- 
ters acting as scaling factors and phase shifts, so called 
extrinsic parameters) can be factorized. This factoriza- 
tion evidences that the signal space can be represented as 
the direct product of two two-dimensional spaces i.e., the 
GW polarization plane and the chirp plane. This repre- 
sentation forms the backbone of the coherent detection 
scheme that follows in the subsequent section. 



In Sec. IV based on a detailed geometrical argument, 
we show that the above signal representation manifests 
the possible degeneracy of the response. This degeneracy 
has been already noticed and studied at length in the 
context of burst detection EH EQ- We investigate 
this question in the specific context of chirps and obtain 
similar results as were presented earlier in the literature. 

In Sec. [V] we obtain the expression of the network 
statistic. Following the principles of the generalized like- 
lihood ratio test (GLRT), the statistic is obtained by 
maximizing the network likelihood ratio over the set of 
unknown parameters. We perform this maximization in 
two steps. We first treat the linear part of the parame- 
terization and show that such a maximization is nothing 
but a least square problem over the extrinsic parameters. 
The solution is obtained by projecting the data onto the 
signal space. We further study the effect of the response 
degeneracy on the resulting parameter estimates. 

The projection onto the signal space is a combined pro- 
jection onto the GW polarization and chirp planes. The 
projection onto the first plane generates two synthetic 
streams which can be viewed as the output of "virtual 
detectors". The network statistic maximized over the ex- 
trinsic parameters can be conveniently expressed in terms 
of the processing of those streams. In practice, the syn- 
thetic streams linearly combine the data from each detec- 
tor in such a way that the GW signature received by each 
detector is added constructively. With this rephrasing, 
the source location angles can be searched over efficiently 

Along with the projection onto the GW polarization 
plane, we also examine the projection onto its comple- 
ment. While synthetic streams concentrate the GW con- 
tents, the so called null streams produced this way com- 
bines the data such that the GW signal is canceled out. 
The null streams are useful to veto false triggers due to 
instrumental artifacts (which don't have this cancellation 
property) . The null streams we obtained here are identi- 
cal to the ones presented earlier in GW burst literature 

In Sec. |VI| we perform the second and final step of 
the maximization of the network statistic over the chirp 
phase function. This step is the difficult part of the prob- 
lem. For the one detector case, we have proposed an ef- 
ficient method, the BCC algorithm which addresses this 
question. We show that this scheme can be adapted to 
the multiple detector case in a straightforward manner, 
hence refer to this as Best Network CC (BNCC). 

Finally, Sec. |VII| presents a proof of principle of the 
proposed method with a full-sky blind search in a sim- 
plified situation. 



II. GENERIC GW CHIRPS 
A. Motivation 

Known observable GW sources e.g., stellar binary sys- 
tems, accreting stellar systems or rotating stars, com- 
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monly involve either orbiting or spinning objects. It is 
not unreasonable to assume that the similar holds true 
even for the unknown sources. 

The GW emission is essentially powered by the source 
dynamics which thus determines the shape of the emit- 
ted waveform. Under linearized gravity and slow mo- 
tion (i.e., the characteristic velocity is smaller than the 
speed of light) approximation, the quadrupole formula 
[TS] predicts that the amplitude of emitted GW is pro- 
portional to the second derivative of the quadrupole mo- 
ment of the physical system. When the dominant part 
of the bulk motion follows an orbital/rotational motion, 
the quadrupole moment varies quasi-periodically, and so 
is the GW. 

The more the information we have about the GW sig- 
nal, the easier/better is the detection of its signature 
in the observations. Ideally, this requires precise knowl- 
edge of the waveform, and consequently requires precise 
knowledge of the dynamics. This is not always possible. 
In general, predicting the dynamics of GW sources in the 
nearly rclativistic regime requires a large amount of ef- 
fort. This task may get further complicated if the under- 
lying source model involves magnetic couplings, mass ac- 
cretion, density-pressure-entropy gradients, anisotropic 
angular momentum distribution. 

Here, we are interested in GW sources where the mo- 
tion is orbital/rotational but the astrophysical dynamics 
is (totally or partially) unknown. While our primary tar- 
get is the unforeseen sources (this is why we remain inten- 
tionally vague on the exact nature of the sources) , several 
identified candidates enter this category because their dy- 
namics is still not fully characterized. These include (see 
[4] for more details and references) binary mergers, quasi 
normal modes from young hot rotating NS, spinning BH 
accreting from an orbiting disk. As motivated before, fol- 
lowing the argument of quadrupole approximation, the 
GW signature for such sources is not completely unde- 
termined: it is expected to be quasi-periodic, possibly 
frequency modulated GW ; in brief, it is a GW chirp. 
This is the basic motivation for introducing a generic 
GW chirp model, as described in the next section. 



B. Generic GW chirp model 

In this section, we describe the salient features of the 
generic GW chirp model used in this paper. We motivate 
the nature of GW polarization, the regularity of its phase 
and frequency evolution. 



1. Relation between the polarizations 

The GW tensor (in transverse traceless (TT) gauge), 
associated to the GW emitted from slow-motion, weak 
gravity sources are mostly due to variations of the mass 
moments (in contrast to current moments) and can be 



expanded in terms of mass multipole moments as [To] 

OO / W 

h TT it) cx £ £ {VVY lm ) STF ^I lm {t - r/c). (1) 
1=2 m=-l 

Here, STF means "symmetric transverse-tracefree" , 
Y lm are the spherical harmonics, and I lm are the mass 
multipole moments. 

We consider here isolated astrophysical systems 
with anisotropic mass distributions (e.g., binaries, ac- 
creting systems, bar/fragmentation instabilities) orbit- 
ing/rotating about a well-defined axis 1 . It can be shown 
that these systems emit GW predominantly in the / = 
|m| = 2 mode. Contributions from any other mass mo- 
ments are negligibly small 2 . The pure-spin tensor har- 
monic (VVy' m ) term provides the GW polarization. 
For / = \m\ = 2 mode, we have 

(VVF 22 ) STF cx (1 + cos 2 e)e+ + 2i cos e e x . (2) 

The tensors e+ and e x form a pair of independent and 
linear-polarization GW tensors (e x is rotated by and an- 
gle 7r/4 with respect to e+). The orbital inclination angle 
e is the angle between the line of sight to the source (in 
Earth's frame) and the angular momentum vector (or the 
rotation axis) of the physical system, see Fig. [l](a). This 
shows that the emitted GW in the considered case carries 
both GW polarizations. 

The GW tensor is fully described as h TT (t) = 
h+(t)e+ + h x (t)e x . The phase shift between the two 
polarizations h + and h x arises from the / term which 
is proportional to the moment of inertia tensor for I = 
\m\ = 2. The quadratic nature of the moment of inertia 
tensor introduces a phase shift of ir/2 between the two 
polarizations h + and h x . This leads to the chirp model 
below 

1 -\- cos 2 6 

h + (t) = A cos(<^(t - t ) + <fo) , (3a) 

h x (t) = Acos e sm(tp(t — to) + 4>o), (3b) 

with to < t < to+T and /i+ jX (t) — outside this interval. 
The phase (f>o is the signal phase at t = to . 

Here, we assume the GW amplitude A to be constant. 
This is clearly an over-simplified case since we indeed 
expect an amplitude modulation for real GW sources. 



1 This condition can be relaxed to precessing systems provided 
that the precession is over time scales much longer than the ob- 
servational time, typically of order of seconds. 

2 Recently, numerical relativity simulations |17| demonstrated that 
this is a fairly robust statement in the specific context of inspi- 
ralling BH binaries. The simulations show that BH binaries emit 
GW dominantly with I = \m\ = 2. However, as the mass-ratio 
decreases, higher multipoles get excited. A similar claim was also 
made in the context of quasi-normal modes produced in the ring- 
down after the merger of two BH, on the basis of a theoretical 
argument, see 118) . page 4538. 
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However we wish here to keep the model simple in or- 
der to focus the discussion on the aspects related to the 
coherent analysis of data from multiple detectors. We 
postpone the study of amplitude modulated GWs to fu- 
ture work. 

The chirp model described in Eq. (3j) clearly depends 
upon several unknown parameters (which need to be es- 
timated from the data) which include the amplitude A, 
the initial phase </>o, the arrival time to of the chirp and 
the inclination angle e. As no precise assumption on the 
exact nature and dynamics of the GW source is made, 
we consider the phase evolution function ip(-) to be an 
unknown parameter of the model ([3| as well. Clearly, it 
is a more complicated parameter than the others which 
are simply scalars. Just like any scalar parameter can be 
constrained to a range of values (e.g., A > 0), the phase 
function ip(-) has to satisfy conditions to be physically 
realistic which we describe in the next section. 



2. Smoothness of the phase evolution 

As explained above, the chirp phase is directly related 
to the orbital phase of the source. The regularity of the 
orbital phase can be constrained by the physical argu- 
ments: the orbital phase and its derivatives are continu- 
ous. The same applies to the chirp phase and derivatives. 

The detectors operate in a frequency window limited 
in the range from few tenths of Hz to a kHz and they 
are essentially blind outside. This restricts our interest 
to sources emitting in this frequency range, which results 
in lower and upper limits on the chirp frequency v(i) = 
(27r) _1 dtp I 'dt and thus on the variations of the phase. 

In addition, the variation of the frequency (the chirp- 
ing rate) can be connected to the rate at which the source 
loses its energy. For isolated systems, this is clearly 
bounded. This argument motivates the following bounds 
on the higher-order derivatives of the phase: 



dp 
~dt 



d 2 v 



dt 2 



< F". 



(4) 



In this sense, Eq.Q determines and strengthens the 
smoothness of the phase/frequency evolution. This is 
the reason why we coined the term "smooth GW chirp" 
in jl]. The choice of the allowed upper bounds F' and 
F" may be based on general astrophysical arguments on 
the GW source of interest. 



III. RESPONSE OF A NETWORK OF 
DETECTORS TO AN INCOMING GW 



A. Coordinate frames 

We follow the conventions of [7] and introduce three 
coordinate frames, namely, the wave frame, the Earth 
frame and the detector frame as given below, see Fig. [T] 



• the wave frame x,. 



(x w ,y w ,z w ) is the frame 



associated to the incoming GW with positive z w - 
direction along the incoming direction and x w — y w 
plane corresponds to the plane of the polarization 
of the wave. 



• the Earth frame x# = (xetVe, ze) is the frame 
attached to the center of the Earth. The xe axis is 
radially pointing outwards from the Earth's center 
and the equatorial point that lies on the meridian 
passing through Greenwich, England. The ze axis 
points radially outwards from the center of Earth 
to the North pole. The xje axis is chosen to form 
a right-handed coordinate system with the x e and 
ze axes. 

• the detector frame = (xd,yd, Zd) is the frame 
attached to the individual detector. The {xd — yd) 
plane contains the detector arms and is assumed 
to be tangent to the surface of the Earth. The Xd 
axis bisects the angle between the detector's arms. 
The Zd axis points towards the local zenith. The 
direction of the yd axis is chosen so that we get a 
right-handed coordinate system. 

A rotation transformation between the coordinate sys- 
tems about the origin is specified by the rotation operator 
O which is characterized by these three Euler angles. We 
define these angles using the "x-convention" (also known 
as z — x — z convention) fTU] . 

Let (0 e , 9 e ,ip e ) and (a, (3, 7) be Euler angles of the ro- 
tation operator relating pairs of the above coordinate sys- 
tems as follows 



x w = 0(cj) e ,9 e ,tp e )x E , 
x d = 0(a,/3,7)x B . 



(5) 
(6) 



All the angles in the Eqs. Q and ^ are related 
to physical/geometrical quantities described in Fig. [T| 
More specifically, we have 



4>e — <j> - 7r/2, 



7T-0, 



= */>, (7) 



After describing the chirp model in the previous sec- 
tion, in this section we derive the response of a network 
of interferometric ground based detectors with arbitrary 
locations and orientations to an incoming GW chirp. The 
first step towards this is to identify the coordinate frames. 



where </> and 9 are the spherical polar coordinates of the 
source in the Earth's frame and the angle vb is the so- 
called polarization-ellipse angle which gives the orienta- 
tion of the source plane. Throughout the paper, we shall 
use 9 and <fi to indicate the source location. 
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FIG. 1: Color online — Coordinate transformations, (a) 0(4> e ,9e,i>e): Earth frame x_b : (xe,Ve , Ze) —* wave frame 
x w : (x m ,y w , z w ), (b) 0(a, f3, 7): Earth frame x_b : (xe, Ve, ze) — * detector frame: x,i : (xd,yd, Zd). The latitude I and longitude 
L of the detector are related to the Euler angles by Eqs. |8| and Q. 



TABLE I: Location and orientation informations of the Earth-based interferometric GW detectors. The location 
of the corner station (vertex) of each detector is given in terms of the latitude and longitude. The orientation of an arm is given 
by the angle through which one must rotate it clockwise (while viewing from top) to point the local North. The corresponding 
detector Euler angles (a, fj, 7) are listed. 



Detector 

TAMA-300 (T) 
GEO-600 (G) 
VIRGO (V) 
LIGO Hanford (LH) 



Vertex 
latitude (N) 

35°40'35.6" 
52°14'42.528" 
43°37'53.0921" 
46°27'18.528" 
LIGO Livingstone (LL) 30°33' 46.4196" 



Vertex 
longitude (E) 

139°32'9.8" 
9°48'25.894" 
10°30'16.1878" 
-119°25'27.5657" 
-90°46'27.2654" 



Arm 1 

90.0° 
25°56'35" 



Arm 2 a 

a-2 

180.0° 229.54° 
291°36'42" 99.81° 
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54.32° 225° 

37.75° 68.778° 

341°25'57.2" 70°34'2.8" 100.5° 46.37° 116.0° 

35.9994° 125.9994° -29.41° 43.55° 171.0° 

107.7165° 197.7165° -0.77° 59.44° 242.72° 



The detector Euler angles (a, j3, 7) are directly related 
to the location and orientation of the detector as follows: 



(8) 
(9) 



a = 


L + n/2, 




13 = 


tt/2 - I, 




1 = 


ai+a 2 




2 


Y 




ai+a 2 


n 




2 


2 



if \ai - o 2 1 > 7T, (10) 
if ja x - a 2 j < 7T, (11) 



where I and L are the latitude and longitude of the corner 
station. The angles a\ and a 2 describe the orientation of 
the first and second arm respectively. It is the angle 
through which one must rotate the arm clockwise (while 
viewing from top) to point the local North. In Table|TJ we 
tabulate the currently running interferometric detectors 
along with their Euler angles. 

Combining Eqs. ^ and we obtain the coordinate 
transformation from the wave frame to the detector frame 
as follows 



where ff e , f e ) = 0(<k, 0e, ^e)0-\a, /?, 7). 



B. Network response 

An interferometric response to the incident GW is ob- 
tained by contracting a GW tensor with the detector 
tensor [see Appendix [B] , which can be re-expressed as 
a linear combination of the two polarizations and h x 
i.e. 



s = f+h+ + f x hy 

= $t[f*h], 



(13) 



0(CM)x d , 



(12) 



and the linear coefficients f + ((j)' e ,0' e ,ijj' e ) and 
f x (<l>' e ,0' e ,ij' e ), commonly termed as the detector 
antenna pattern functions, represent the detector's 
directional response to the + and x polarizations 
respectively. For the compact expression provided by 
Eq. (13 1, we have defined the complex GW signal to be 
h = hj- + ih x and the complex antenna pattern function 
f++ifx- 



h + 
to be / 
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The detector response s and the incident GW signal 
h are both times series. In a network where the various 
detectors are located at different locations on the Earth, 
for example the LIGO- Virgo network, the GW will arrive 
at the detector sites at different time instances. However, 
all the measurements at the various detectors need to be 
carried out with a reference time. Here, just for our con- 
venience, we choose the observer attached to the Earth's 
center as a reference and the time measured according to 
this observer is treated as a reference. Any other refer- 
ence would be equally acceptable. 

Assuming that the detector response is labeled with 
this time reference (a single reference time for all detec- 
tors in a detector network), we have 

S (t) = viirh(t-T(^6))}, (u) 

where t(</), 9) = (r^— r£)-w(</>, 9)/c denotes the difference 
in the arrival times of the GW (propagating with the unit 
wave vector w) at the detector and at the center of the 
Earth located at and rg, respectively. Note that this 
value can be positive or negative depending on the source 
location. 



C. Vector formalism 

In the following, we distinguish scalars by using roman 
letters, vectors are denoted by small bold letters, and 
matrices by bold capitals. We denote the fc-th element 
of vector a by a[k] and correspondingly, the element of 
matrix A at row k and column I by A[k, I]. The matrices 
A T and A H = (A T )* designate the real and hermitian 
transposes of A respectively. 

We consider now a GW detector network with d inter- 
ferometers. Each detector and its associated quantities 
are labeled with an index j = 1, . . . , d which we also use 
as a subscript if required. We assume that the output 
response of each detector is sampled at the Nyquist rate 
v s = l/t s where t s is the sampling interval. We then 
divide the data in blocks of N consecutive samples. In 
this set-up, the detector as well as the network response 
is then defined by forming vectors with these blocks of 
data. 

Let us consider a given GW chirp source at sky location 
(<j), 9) . Let the response of the j-th detector be Sj with 
entry Sj[k] — Sj{tk + Tj((f>,9)), where tk — to + kt s , k — 
0, . . . , N — 1, and t is the reference time i.e. the time of 
arrival of GW at the center of the Earth. Note that, with 
the above definition, we compensate for the time delay 
Tj((f>,9) between the detector j and the Earth's center. 
Thus, in this set-up, the GW signal starts and ends in 
the same rows in the data vectors Sj of all the detectors. 

For compactness, we stack the data from all the detec- 
tors in the network into a single vector s of size Nd x 1 , 
such that s T = [s^sj . . . sj] forms the network response. 
In this convention, the network response can be expressed 
compactly as the Kronecker product (see Appendix |A| for 



the definition) of the network complex beam pattern vec- 
tor f = {fj, j = 1 . . .d} G C dxl and the complex GW 
vector h = {h(t k ),t k = t + t s k with k = . . . N - 1} G 
C Nxl viz, 

s = »[f*®h]. (15) 

The above expression is general enough to hold true for 
any type of incoming GW signal. The Kronecker prod- 
uct in this expression is the direct manifestation of the 
fact that the detector response is nothing but the tensor 
product between the detector and the wave tensors. 

D. GW chirp as a linear model of the extrinsic 
parameters 

In previous section, we have obtained the network re- 
sponse to any type of incoming GW with two polariza- 
tions. In what follows, we wish to investigate how this re- 
sponse manifests in case of a specific type of GW, namely 
GW chirp described in Eq. pi. We also want to under- 
stand how various parameters explicitly appear in the 
network response. 

It is insightful to distinguish the signal parameters 
based on their effect on the signal model. The param- 
eters are separated into two distinct types traditionally 
referred to as the "intrinsic" and "extrinsic" parameters. 
The extrinsic parameters are those that introduce scaling 
factors or phase shifts but do not affect the shape of the 
signal model. Instead, intrinsic parameters significantly 
alter the shape of the signal and hence the underlying 
geometry. 

The network response s mingles these two types of 
parameters. Our work is considerably simplified if we 
can "factorize" the extrinsic parameters from the rest. 
For the chirp model described in Eq. pj, we count four 
extrinsic parameters, namely {A, <pQ, e, ip e } and perform 
this factorization in two steps. 

1. Extended antenna pattern includes the inclination angle 

We absorb the inclination angle e into the antenna pat- 
tern functions and rewrite the network signal as 

s = R|f*<3h], (16) 

where h = ae is the GW vector. It only depends on the 
complex amplitude a = A exp i(f>o and on the phase vector 
e = {cxp (i (p[k]) , with ip[k] = <p(kt s ), k = . . . N — 1}. 

The extended antenna pattern f incorporates the incli- 
nation angle e as follows 3 

- 1 + cos 2 e 

1 = r + + i cos e r x . (17) 



We remind the reader that similar quantity was previously in- 
troduced in Eq. (3.19) of [7]. 
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2. Gel' f and functions factorize the polarization angles from 
source location angles 

The second step is to separate the dependency of f on 
the polarization angles {^i 6 } from the source location 
angle and the detector orientation angles. The earlier 
work 7 shows that the Gel'fand functions (which are 
representation of the rotation group 50(3)) provide an 
efficient tool to do the same. For the sake of complete- 
ness, Appendix [B] reproduces some of the calculations of 
0. The final result (see also Eqs. (3.14-3.16) of [7 ) yields 
the following decomposition: 



columns of E which we shall refer to as chirp plane. These 
two spaces embody two fundamental characteristics of 
the signal: the former characterizes gravitational waves 
while the latter characterizes chirping signals. The Kro- 
necker product in the expression of II shows explicitly 
that the network response is the result of the projection 
of incoming GW onto the detector network. 

The norm of the network signal gives the "signal" (and 
not physical) energy delivered to the network, which is 



f = t+d + t-d* 



(18) 



where the vector d £ C dxl carries the information of the 
source location angles (<f>,0) via (<j) e ,9 e ) and the detector 
Eulcr angles {a 3 ■, (3j, jj}. Its components are expressed 
as 



d\j] 



n=-2 



The tensor T mn designates rank-2 Gel'fand functions. 
The coefficients t+ and t_ depend only on the polar- 
ization angles {ip,e}, viz. 

t± = T 2±2 (V, e, 0) = ^ — -L exp( T 2iV)- (20) 



Finally, we combine Eqs. (16 1 and (|18|) and obtain an 



expression of the network response where the extrinsic 
parameters are "factorized" as follows, 



n P . (2i) 






( at*_\ 




a*t + 




a t* + 




{ a*t_ J 


p 



Eq. (21 1 evidences the underlying linearity of the GW 
model with respect to the extrinsic parameters. The 4- 
dimensional complex vector p defines a one-to-one (non- 
linear) mapping between its components and the four 
physical extrinsic parameters {A, 4>q, e, i/j} (we will detail 
this point later in Sec. VA3). Note that the first and 



fourth components as well as the second and third com- 
ponents of p are complex conjugates. This symmetry 
comes from the fact that the data is real. 

The signal space as defined by the network response 
is the range of II and results from the Kronecker prod- 
uct of two linear spaces: the plane of C d generated by 
the columns of D which we shall refer to as GW po- 
larization plane 4 and the plane of C N generated by the 



4 In [7] , this plane was referred to as "helicity plane" because it is 
formed by the network beam patterns for all possible polariza- 
tions. 



NA 2 



(22) 



, d e ,0)[T 2n (aj, /3j7j)-T_2n(aj, Pj,lj)Y 

(19) 



Clearly, the dependence on the number of samples N 
implies that the longer the signal duration, the larger is 
the signal energy and is proportional to the length of the 
signal duration. The factor ||f || is the modulus of the 
extended antenna pattern vector. It can be interpreted 
as the gain or attenuation depending on the direction of 
the source and on the polarization of the wave. 



IV. INTERPRETATION OF THE NETWORK 
RESPONSE 



In this section, we focus on understanding the underly- 
ing geometry of the signal model described in Eq. (21 1. A 



useful tool to do so is the Singular Value Decomposition 
(SVD) [20 . It provides an insight on the geometry by 
identifying the principal directions of linear transforms. 



A. Principal directions of the signal space: 
Singular Value Decomposition 

The SVD is a generalization of the eigen-decomposition 
for non-square matrices. The SVD factorizes a matrix 
A e C mxn into a product A = \] A H A W^ of three ma- 
trices U A e C mw , S A £ R rxr and V A e C nxr where 
r < to, n is the rank of A. The columns of and V^ 
are orthonormal i.e., U^U^ = V^Va = I r . The diago- 
nal of Ha are the singular values (SV) of A. We use here 
the so-called "compact" SVD (we retain the non-zero SV 
only in the decomposition), such that the matrix Ha is 
a positive definite diagonal matrix. 

The SVD is compatible with the Kronecker product 
[21] : the SVD of a Kronecker product is the Kronecker 
product of the SVDs. Applying this property to II, we 
get 



n= (U D ®U E )CZ D ®-Z E )(V D ®V E ) H . (23) 

Therefore, the SVD of II can be easily deduced from 
the one of D and E. We note that D and E have simi- 
lar structure (two complex conjugated columns), see Eq. 
(21 1. In Appendix [C] we analytically obtain the SVD 
of a matrix with such a structure. Thus, applying this 
result, we can straightaway write down the SVDs for D 
and E as shown in the following sections. 
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1. GW polarization plane: SVD of D 
Let us first introduce some variables 

d 

V ee & H & = Y J \d[3\\\ 

i=i 

A EE d T d = Y / d[j)\ 

3=1 

(5 ee argA. 

In the nominal case, the matrix D has rank 2, viz. 



J D 



o i 

CT 2 



(24) 

(25) 
(26) 
i. 

(27) 



with two non-zero SV o\ = \JT> + | A| and 02 = 
y/V — |A| ((Ti > cr 2 ) associated to a pair of left-singular 
vectors = [vi,V2] with 



vi 



exp(— id) 
1 



v 2 



exp(— iS) 
-1 



and of right-singular vectors Uu = [111,112] with 



ui 



u 2 



exp(— id)d + d* 
exp(— iS)d — d* 



(28) 

(29) 
(30) 



Note that the vector pair {ui,u 2 } results from the 
Gram-Schmidt ortho-normalization of {d, d*}. 

Barring the nominal case, for a typical network built 
with the existing detectors and for certain sky locations 
of the source, it is however possible for the smallest SV 02 
to vanish. In such situation, the rank of D reduces to 1. 
We then have T, D = cr 1: V D = vi and = m. We give 
an interpretation of this degeneracy later in Sec. |IV B| 



2. Chirp plane: SVD of E 

The results of the previous section essentially apply to 
SVD calculation of E. However, there is an additional 
simplification due to the nature of the columns of E. 
Indeed, the cross-product 



the case of interest. We can thus consider 5 that e T e « 
and e H e = N. Therefore, following Appendix [c| 2, the 
SVD of E is given by S B = y/Nl 2 /2, V E = I 2 'and 
U s = 2E/VN. 



3. Signal space: SVD of II 
We obtain the SVD for II using the compatibility of 



the SVD with the Kronecker product stated in Eq. ( 23 1 
In the nominal case where D has rank 2, we have 



N 
2~ 



(7ll 2 2 
O2 O2I2 



with four left-singular vectors 
Vi <g) I 2 

and four right-singular vectors 
Un 



V2 



(32) 



(33) 



2 

N 



Ui e Ui 



u 2 (g) e u 2 



(34) 



B. The signal model can be ill-posed 

In the previous section, we obtained the SVD of II in 
the nominal case where the matrix D has 2 non-zero SVs. 
As we have already mentioned, for a typical detector net- 
work, there might exist certain sky locations where the 
second SV o~2 of D vanishes which implies that the rank 
of D degenerates to 1. In such cases, this degeneracy 
propagates to II and subsequently its rank reduces from 
4 to 2. 

In order to realize the consequences of this degeneracy, 
we first consider a network of ideal GW detectors (with 
no instrumental noise). Let a GW chirp pass through 
such a network from a source in a sky location where 02 = 
0. The detector output is exactly equal to s. An estimate 
of the source parameters would then be obtained from 
the network data by inverting Eq. (21 1. However, in this 
case, this is impossible since it requires the inversion of 
an under-determined linear system (there are 4 unknowns 
and only 2 equations). 

This problem is identical to the one identified and dis- 
cussed at length in a series of articles devoted to un- 
modeled GW bursts [3 HH1 HI], where this problem is 
formulated as follows: at those sky locations where D 



N-l 

e T e = ^2 exp(2iy[fc]) , (31) 

k=0 

is an oscillating sum. This sum can be shown [3] to be of 
small amplitude under mild conditions compatible with 



5 This amounts to saying that the two GW polarizations (i.e., the 
real and imaginary parts of expiip[k]) are orthogonal and of equal 
norm. Note that this approximation is not required and can be 
relaxed. This would lead to use a version of the polarization pair 
ortho-normalizcd with a Gram-Schmidt procedure. 
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is degenerated, the GW response is essentially made of 
only one linear combination of the two GW polariza- 
tions. It is thus impossible to separate the two individual 
polarizations (unless additional prior information is pro- 
vided). We want to stress here that this problem is not 
restricted to unmodeled GW bursts but also affects the 
case of chirping signals (and extends to the chirps from 
inspiralling binaries of NS or BH 6 ). This is mainly be- 
cause the degeneracy arises from the geometry of the GW 
polarization plane which is same for any type of source. 

The degeneracy disappears at locations where 02 > 
even if it is infinitcsimally small. However, when 02 is 
small, the inversion of the linear equations in Eq. (21 1 is 
very sensitive to perturbations. With real world GW de- 
tectors, instrumental noise affects the detector response 
i.e., perturbs the left-hand side of Eq. 



(a) LH + LL 



(b)LH + LL + V 



(211 



A useful tool to investigate this is the condition number 
llj. It is a well-known measure of the sensitivity of linear 
systems. The condition number of a matrix A is defined 
as the ratio of its largest SV to the smallest. For unitary 
matrices, cond (A) = 1. On the contrary, if A is rank 
deficient, cond (A) — > 00. For the matrix II, we have 



cond (II) 



cr 2 



lv + 


|A| 


V- 


|A| 



(35) 



In Fig. [2] we show full-sky plots of 1 /cond (II) for vari- 
ous configurations (these figures essentially reproduce the 
ones of [H]). We see that, even for networks of mis- 
aligned detectors, there are significantly large patches 
where cond (II) takes large values. In those regions, the 
inversion of Eq. (21 1 is sensitive to the presence of noise 



and the estimate of the extrinsic parameters thus have a 
large variance. 

a. Connection to the antenna pattern function 
Interestingly, the SVs of D and F = [f f *] coincide. This 
can be seen from the following relationships we directly 
obtained from the definitions in Eqs. (17 1, (18) and (20) 



f F 



D 



t + t*_ 




"l*+l \t-\ 




= F 




t- t* + 




\t-\ IM 



(36) 



f 2 

9 1 



Iff 



0.8 
0.6 
0.4 

0.2 



(radian) 
(c) LH + LL + V + G 



ij> (radian) 
(d) LH + LL + V + G + T 



-202 

(radian) 



-2 2 

i|> (radian) 



FIG. 2: Degeneracy of the network response. We show 
here the inverse of cond (II) for various detector networks (the 
abbreviated detector names are listed in Table[I]). The brighter 
regions of the sky correspond to the large conditioning number 
cond (II). The fraction of the sky where 1/cond (II) < 0.1 is 
(a) 25% (b) 4.5% (c) ~ 0% (d) ~ 0%. Since the LIGO detec- 
tors are almost aligned and they show the largest percentage 
of degeneracy. 



When a 2 ~ 0, we thus have |f+ x f x | 2 <~ where 
f + = 3?[f] and f x = 3[f] are the network antenna pat- 
tern vectors. This means that at such sky locations, the 
antenna pattern vectors get aligned even if the detectors 
in the network are misaligned. In other words, despite 
the considered network is composed of misaligned detec- 
tors, it acts as a network of aligned detectors at those sky 
locations. (Of course, for perfectly co-aligned detectors, 
f + cx f x at all sky locations.) Networks with many differ- 
ent detectors having different orientations are less likely 
to be degenerate. This is confirmed in Fig. [2] where we 
see that the size of the degenerate sky patches reduces 
when the number of detector with varied orientations in- 
creases. For instance, with a network of LL-LH-G-V, 
(assuming they have identical noise spectrum), it goes to 
zero. 



The matrix F can be obtained from D by a unitary 
transformation. Both matrices share the same singular 
spectrum. We can therefore write 



- 2 2 = £l/bH 2 ~ 



£/b1 



(37) 



6 Contrary to the generic chirp model considered here, the phase 
and amplitude functions of inspiralling binary chirps follow a pre- 
scribed power-law time evolution. These differences affect only 
the geometry of the "chirp plane", but not that of the "GW 
polarization plane", hence the conclusion on the degeneracy re- 
mains the same. 



V. NETWORK LIKELIHOOD ANALYSIS: GW 
POLARIZATION PLANE AND SYNTHETIC 
STREAMS 

Generally speaking, a signal detection problem 
amounts to testing the null hypothesis (Hq) (absence 
of signal in the data) vs the alternate hypothesis (Hi) 
(presence of signal in the data). Due to the presence 
of noise, two types of errors occur: false dismissals (de- 
cide Hq when Hi is present) and false alarms (decide 
Hi when Hq is present). There exist several objective 
criteria to determine the detection procedure (or statis- 
tic) which optimizes the occurrence of these errors. We 
choose the Neyman-Pearson (NP) approach which mini- 
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mizes the number of false dismissals for a fixed false alarm 
rate. It is easily shown that for simple problems, the like- 
lihood ratio (LR) is NP optimal. However, when the sig- 
nal depends upon unknown parameters, the NP optimal 
(uniformly over all allowed parameter values) statistic is 
not easy to obtain. Indeed for most real- world problems, 
it does not even exist. However, the generalized likeli- 
hood ratio test (GLRT) [22] have shown to give sensible 
results and hence is widely used. In the GLRT approach, 
the parameters are replaced by their maximum likelihood 
estimates. In other words, the GLRT approach uses the 
maximum likelihood ratio as the statistic. Here, we opt 
for such a solution. 

As a first step, we consider the simplified situation 
where all detectors have independent and identical in- 
strumental noises and this noise is white and Gaussian 
with unit variance. We will address the colored noise case 
later in Sec. EE 

In this case, the logarithm of the network likelihood 
ratio (LLR) is given by 



A(x) = -llx - s] 



(38) 



where || • || is the Euclidean norm (here in R Nd ) and we 
omitted an unimportant factor 1/2. The network data 
vector x is constructed on the similar lines as that of the 
network response s, i.e. first, it stacks the data from 
all the detectors into x T = [xf , xj, . . . xj] and then at 
each detector, the data is time-shifted to account for the 
delay in the arrival time Xj = {x-Jfc] = xj(tk + Tj),t^ — 
t s k and k = 0...N -1}. 



A. Maximization over extrinsic parameters: scaling 
factors and phase shifts 

Following the GLRT approach, we maximize the net- 
work LLR A with respect to the parameters of s. We 
replace s by its model as given in Eq. (21 1 and consider 



at first the maximization with respect to the extrinsic 
parameters p. 

1. Least-square fit 

The maximization of the network LLR over p amounts 
to fitting a linear signal model to the data in least square 
(LS) sense, viz. 

I|2 II „ ||2 



minimize 



A(x) 



x — lip over p. (39) 



This LS problem is easily solved using the pseudo- 
inverse II# of II [20] . The estimate of p is then given 

by 

p = n # x. (40) 

The pseudo-inverse can be expressed using the SVD 
of II as II# = VnSn^Urf (note that II# is always de- 
fined since we use the compact SVD restricted to non-zero 
SVs). 



Substituting Eq. (40 1 in Eq. ( 39 ) , we get the LS mini- 
mum to be 



A(x) 



|x-U n U£x|| 2 , 



(41) 



where we used V^Vn = I r . Eq. (41 1 can be further 
simplified into 7 



A(x) 



|Uf[x| 



(42) 



It is interesting to note that the operator UnU^ is a 
(orthogonal) projection operator onto the signal space 
(over the range of II) i.e. U n U^II = (nn#)n = II. 



2. Signal-to-noise ratio 

The signal-to-noise ratio (SNR) measures the level of 
difficulty for detecting a signal in the noise. In the present 
case, along with the amplitude and duration of the in- 
coming GW, the network SNR also depends on the rel- 
ative position, orientation of the source with respect to 
the network. Therefore, the SNR should incorporate all 
these aspects. A systematic way to define the SNR is to 
start from the statistic. 

Let the SNR p of an injected GW chirp Sq = IIpo be 



P 2 ^A(s ). 



(43) 



Note that in this expression, the matrix II in the statis- 
tic and in So are the same. Using the SVD of II and the 
property of the projection operator U^Urr = Ir, we get 

p 2 = || So ||". The SNR is equal to the "signal energy" in 
the network data as defined in Eq. (22 1. Thus, the SNR p 
scales as 



'N as expected and it depends on the source di- 
rection, polarization and network configuration through 
the gain factor ||f|| 9 . Fig. 3 illustrates how this factor 
varies for the network formeaby the two LIGO detectors 
and Virgo. Fig. [3] displays the ratio p/pbest between the 
global SNR (obtained with a coherent analysis) and the 
largest individual SNR (obtained with the best detector 
of the network). The panels (a) and (b) are associated 
to the "worst" (minimum over all polarizations angles e 
and ip) and "best" (maximum) cases respectively. Ide- 
ally, when the detectors are aligned, the enhancement 
factor is expected to be \fd (w 1.73 in the present case). 
In the best case, the enhancement is > 1.7 for more than 
half of the sky (94% of the sky when > 1.4). In the worst 
case, the SNR enhancement is 1.28 at most and 8.5% of 
the sky gets a value > 1.1. 



7 For inspiral case, this expression is equivalent to Eq. (4.8) of [JJ. 

8 If the noise power is not unity, it would divide the signal energy 
in this expression. When we have only one detector, the SNR p 2 
is consistent with the definition usually adopted in this case. 

9 The SNR p 2 is similar to b 2 defined in Eq. (3.17) of [7] in case 
of inspiralling binary signal and colored noise. 
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(a) LH + LL + V: min[p/Best(p.)] 




(b) LH -t LL + V:max|p/Best(p.)] 



1.05 o.S| 



Q> (radian) 



o (radian) 



II: 



FIG. 3: Benefits of a coherent network analysis (SNR 
enhancement). We display the polar maps of the following 
quantities for the LL-LH-V network (a) min^, i(! p/pbcst and 
(b) max^,j p/pbost- Here, pbost designates the best SNR of 
the detectors in the network. The maximum, minimum are 
taken over all the polarization angles {4>, e}. 



3. From geometrical to physical parameter estimates 

The components of p do not have a direct physi- 
cal interpretation but as mentioned earlier, they are 
rather functions of the physical parameters. Following 
the above discussion, if we assume that we obtained pa- 
rameter estimates p from the data through Eq. (40 1, 
then one can retrieve the physical parameters A, (f>o, 



e and tp by inverting the non-linear map which links 

p = (a t*_ a*t + a t* + a*<_) to these parameters as 
given below 



that the estimation goes worst with the conditioning of 

n. 

Regularization techniques seem to give promising re- 
sults in the context of GW burst detection [DJ EH EQ • 
Following this idea, we may consider to "regularize" the 
LS problem in Eq. ( 39 1 . To do so, additional information 



on the expected parameters is required to counterbal- 
ance the rank-deficiency. Unfortunately, we don't expect 
p to follow a specific structure. The only sensible prior 
that can be assumed without reducing the generality of 
the search is that ||p|| is likely to be bound (since the 
GW have a limited amplitude A). It is known [23] that 
this type of prior is associated to the use of the so-called 
Tykhonov regulator and that we don't expect significant 
improvements upon the non-regularized solution. 

Probably, one difference might explain why regulariza- 
tion techniques do not work in the present case while it 
does work for burst detection. We recall that in the burst 
case, the parameter vector comparable to p are the sam- 
ples of the waveform. This vector being a time-series, 
it is expected to have some structure, in particular it is 
expected to have some degree of smoothness. The use of 
this a priori information improves significantly the final 
estimation. 

While regularization will not help for the estimation 
of the extrinsic parameters, they may be of use to im- 
prove the detection statistic. We consider this separate 
question later in Sec. 



VII B 



A = 



- [arg(p[l]) 



arg(p[2])] 



1 



cos 



>rg(p[l])+arg(p[2])] , 



(44) 
(45) 

(46) 
(47) 



4- Degeneracy and sensitivity of the estimate to noise 

Upper-bounds for the estimation error can be ob- 
tained using a perturbative analysis of the LS problem 
in Eq. (|39|. A direct use of the result of [2D], Sec. 5.3.8 
yields 



(48) 



VN 

< - — cond (n). 
P 



This bound is a worst-case estimate obtained when the 
noise term which affects the data x is essentially concen- 
trated along the directions associated to the smallest SV 
of II. The noise is random and it spans isotropically all 
Nd dimensions of the signal space. As described above, 
the space associated to the smallest SV has only 2 di- 
mensions. Therefore, the worst case is very unlikely to 
occur and the above bound is largely over-estimated on 
the average. However, it gives a general trend and shows 



B. Implementation with synthetic streams 

In the previous section, we maximize the network LLR 
with respect to the ex trinsic parameters resulting in the 
statistic A in Eq. (42 1. Here, we obtain a more simple 



and practical expression for A which will be useful for 
maximization over the remaining intrinsic parameters. 
From Eqs. (p3l and d42|, we have 



k(x)=\\(\J D ®XJe) H x\ 



(49) 



It is useful to reshape the network data x into a N x d 
matrix X = [xiX2 . . . x<j]. This operation is inverse to 
the stack operator vec() defined in Appendix |A| 



Using a property of the Kronecker product in Eq. ( A2 1 
we obtain the reformulation 



A(x) = ||vec(UfXU^)| 



(50) 



There are two possibilities to make this matrix prod- 
uct, each being associated to a different numerical imple- 
mentation for the evaluation of A. 

A first choice is to first multiply X by U|( and then by 
UJj. In practice, this means that we first compute the 
correlation of the data with a chirp template, then we 
combine the result using weights (related to the antenna 
pattern functions). This is the implementation proposed 
in [TJ. It is probably the best for cases (like, searches 
of inspiralling binary chirps) where the number of chirp 
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templates is large (i.e., larger than the number of source 
locations) and where the correlations with templates are 
computed once and stored. 

The second choice is to first multiply X by and 
then by which we adopt here. This means that 
we first compute Y = XUJ, which transforms the net- 
work data into two TV-dimensional complex data vectors 
[yi > Yz] = Y through an "instantaneous" linear combi- 
nation. Then, we correlate these vectors with the chirp 
template. We can consider yi and y2 as the output of 
two "virtual" detectors. For this reason, we refer to those 
as synthetic streams in connection to [24 who first coined 
the term for such combinations. Note that, irrespective 
of the number of detectors, one always gets at most two 
synthetic streams. We note that though the synthetic 
streams defined in |24j are ad-hoc (i.e., they have no re- 
lation with the LR) , the ones obtained here directly arise 
from the maximization of the network LLR. 

We express the network LLR statistic in terms of the 
two synthetic streams as 

A(x) = i (|e* yi | 2 + |e T yi | 2 + |e"y 2 | 2 + |e T y 2 | 2 ) , 

(51) 

where y; = Xu ; * for I = 1, 2. This expression can be fur- 
ther simplified by using the symmetry (easily seen from 
Eqs. p9l) and @), 



= exp(i<5)ui 
We finally obtain 



-exp(z<5)u 2 . (52) 



A(x)^(|e" yi | 2 



|e"y 2 | 2 ). 



(53) 



The linear combination in each stream is such that 
the signal contributions from each detector add up con- 
structively. In this sense, synthetic streams are similar 
to beam-formers used in array signal processing |25j . 

When the data is a noise free GW chirp, i.e., x = s, 
we then have 

yi[k] = P T (D <g) E[fc]) T u ; * = p T (a iV; * <g E[fc] T ) . (54) 

Here, E[fc] represents the fc-th row of E. Writing ex- 
plicitly, we have 



Y2 



ia 2 



3feh}e i5/2 ■ 



(55) 



( t _ e -M/2 _ 
(l6l. This 



where q x = {t^ lA ' 2 + t + e ld ' 2 )*, q 2 
t + e l5 / 2 )* and where h is as defined in Eq. 
shows that the synthetic streams y; are rescaled and 
phase shifted copies of the initial GW chirp h which en- 
hances the amplitude of the signal by the appropriate 
factor as shown above. 

To assess this enhancement, we compute the SNR per 
synthetic stream as below. 




(b) LH + LL + V: max[p/p t J 



- 




t 




-2 


2 




(ft (radian) 



FIG. 4: Color online — SNR per synthetic streams 
and benefits of a coherent network analysis (SNR en- 
hancement). We display the polar maps of the following 
quantities for the LL-LH-V network: (a) {pij p\)^^ an d (b) 
maxy, |f pi max; pi . W e denote pi the SNR of synthetic stream 
as defined in Eq. (571. The maximum and average are taken 



over all the polarization angles {ip,e}. 



1. SNR per synthetic streams 



The network SNR can be split into the contribu- 
tions from each synthetic stream i.e. we write p 2 = 



|£ n Vjfp || , as 



2 2 
P = Pi 



pI 



(56) 

where we define pi = VNui || (v^ (g) I2) "^Po || /2 for Z = 1,2. 
More explicitly, we have 



Pi 



<Ji\qi\A. 



(57) 



The synthetic streams contributes differently depend- 
ing on the polarization of the incoming wave. Fig. [^il- 
lustrates this with the network formed by the two LIGO 
detectors and Virgo. 

Let us assume that po is randomly oriented. Since 
Vi and v 2 have unit norms, we get the average value 
< (p2/pi)e,ip 1/cond (II) < 1 for most of the sky 
as indicated in Fig. [5] (a). Note that this panel matches 
well with Fig. [2] (b). Thus, on average, yi contributes 
more to the SNR than y 2 . However, the situation may 
be different depending on the specific polarization state 
of the wave. Fig. [4] (b) show the maximum of the ratio 
pj max; =12 p\ for all polarization angles e and tp mostly 
takes the value » %/2- This implies that for a given di- 
rection, one can always find polarization angles where 
the two synthetic streams contribute equally. This holds 
true for all the sky locations, except at the degenerate 
sky locations where y 2 does not contribute, hence the 
SNR ratio is 1. 

Another way to see this is to examine the expression 
of the SNR difference in terms of the signal and network 
quantities, namely 



Pi 



p\ = 



NA 2 



1 



2 \ 2 
cos e 



f+ - fx 



cos 2 e 



2IAI 



(58) 
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where we have 10 |A| 2 



-) 2 -4||f+xfx 



In the simple case of a face-on source, i.e. e = 0, 
and if f x is orthogonal to f + with 1 1 f x 1 1 = 1 1 f+ L then 
both synthetic streams carry the same SNR p\ — P2 — 
•y/7V/2A||f + ||. However, for other situations, some other 
wave polarization would lead the same. 

Inversely looking at min £i ^ p/max; pi, it can be shown 
that one can always find a GW polarization such that 
one of the synthetic streams does not contribute to the 
SNR. 



C. Null streams 

1. Review and relation to synthetic streams 

The access to noise only data is crucial in signal detec- 
tion problem. Such data is not directly available in GW 
experiments, but the use of multiple detectors allows to 
access it indirectly using the null streams. The general 
idea behind the null stream is to construct a data stream 
from the individual detector streams which nullifies the 
signature of any incoming GW from a particular direc- 
tion. Since this signal cancellation is specific to GWs, 
null streams naturally provide an extra tool to verify that 
a detected signal is indeed a GW or instead a GW like 
features mimicked by the detector noise whose detection 
thus has to be vetoed. This is a powerful check since it 
does not require detailed information about the potential 
GW signal under test, except an estimate of its source 
location. (Note that in practice, the implementation of 
the veto test may be complicated by the imprecision of 
the direction of arrival and of the errors of calibration 
V The existence of null streams has been first iden- 
tified in [T2] in the case of three detector networks. At 
present, handful of literature 0111] exists on the use of 
null streams in GW data analysis. 

Null streams are usually introduced as a general post- 
processing of the data independent of the detection of 
specific GWs. Below, we make this connection in the 
domain of our formalism. We recall that the network 
data at a given time (e.g. , the first row of the matrix 



X introduced in Sec. VBl is a d-dimensional vector in 



K . This space is a direct sum of the GW polarization 
plane and its orthogonal complementary space. We have 
shown that the GW polarization plane is a 2-dimensional 
space, spanned by a pair of orthonormal basis vectors 
which are associated to the two synthetic streams. The 



10 The synthetic streams (on the average sense) are also connected 
to the directional streams introduced in the context of LISA 26 . 
If we integrate pf over the inclination and polarization angles 

€,i>, we obtain <(pf - p 2 2 )/2), ti , = 2|A|/5 and (||f || 2 } E ,,/. = 223/5. 
Thus, the SNRs of the synthetic streams p 2 when averaged over 
the polarization angles are proportional to the SNRs obtained by 
v+ and v x — the directional streams in the LISA data analysis, 
see Eqs. (25-28) of [26]. 



complementary space to the GW polarization plane is 
ad-2 dimensional space and it is spanned by d — 2 
"null vectors" . Similarly to the synthetic streams, the 
null streams can be constructed from these null vectors. 
Thus, the numbers of synthetic and null streams sum 
up to d. Nominally, we have d — 2 null streams. How- 
ever, when the GW polarization plane degenerates to a 
1-dimensional space {a 2 ~ 0) as explained in Sec. IV B 
the number of null vectors becomes d — 1. For a two 
detector network, in the nominal case, there is no null 
stream as d — 2 = 0. However, for degenerate directions, 
one can construct a null stream. For aligned pair of de- 
tectors (as it is almost the case for the two LIGO LH and 
LL), the fraction of the degenerate sky location is large, 
see Fig. [2] This null stream would turn out to be useful 
for vetoing in this case. 

In the next section, we explain how the null streams 
can be obtained numerically in the nominal case. The 
extension to the degenerate case is straightforward. 



2. Obtaining the null streams numerically 

The numerical construction of the null streams can be 
achieved in various ways. One such approach could be to 
obtain the full SVD of D and construct the null streams 
from the eigen-vectors corresponding to the zero SVs. 
This approach was taken in [2]. Here, we take an al- 
ternative approach. We construct the null streams by 
successive construction of orthonormal vectors via multi- 
dimensional cross product as described below. 

Assuming some direction of arrival, we express any in- 
stantaneous linear combination of the time-shifted data 
(to compensate for different time of arrivals at the detec- 
tors' site with respect to the reference) as 



y(x) = Xu, 



(59) 



where the vector u g 



-■dxl 



contains the tap coefficients. 



Eq. ( 59 1 can be rewritten as 

y(x) = vec(lArXu) = (u T <g> I N )x. (60) 

The vector u defines a null-stream if y(x) = On when- 
ever x is a GW. Let us assume that we indeed observe a 
GW chirp i.e., x — Sq = IIpo. We thus have 



y(s ) = [(u T U D ) ®U £ ]£ n V£p . 



(61) 



If u is in the null space of U r> , the null-stream condi- 
tion is satisfied for all po. Since the null space of Ud is 
orthogonal to its range, an obvious choice for u is 



U = Ui X U2 



d* x d 

VT> 2 - A 2 ' 



(62) 



Nominally, Urj is a 2-dimcnsional plane in C d . Its null 
space is therefore d — 2 dimensional. An orthonormal 
basis of this space can be obtained recursively starting 
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from U3 = u as defined above and applying the following 
generalized vector cross-product formula for n > 3: 



u n [i] = eijki... m ui[j)u 2 [k]u 3 [l) . . .u„_i[m]. 



(63) 



Here, . m is the Levi-Civita symbol 11 . The u n de- 
notes an orthonormal set of d — 2 vectors, {u„, for 3 < 
n < d}. The components of these vectors are the tap co- 
efficients to compute the null-streams. By construction, 
the resulting null streams are uncorrclated and have the 
same variance. 

To summarize the main features of our formalism. The 
representation of a GW network response of unmodeled 
chirp as a Kronecker product between the GW polariza- 
tion plane and the chirp plane forms the main ingredi- 
ent of this formalism. Such a representation allows the 
signal to reveal the degeneracy in a natural manner in 
the network response. It also evidences the two facets 
of the coherent network detection problem, namely, the 
network signal detection via synthetic streams and veto- 
ing via null streams. The coherent formalism developed 
in [7] for inspiralling binaries lacked this vetoing feature 
due to the difference in the signal representation. 

In the rest of this paper, we do not dis- 
cuss/demonstrate the null streams applied as a vetoing 
tool to the simulated data. This will be demonstrated in 
the subsequent work with the real data from the ongoing 
GW experiments. 



D. Colored noise 

The formalism developed till now was exclusively tar- 
geted for the white noise case. We assumed that the noise 
at each detector is white Gaussian. In this subsection, 
we extend our formalism to the colored noise case. We 
remind the reader that the main focus of this paper is to 
develop the coherent network strategy to detect unmod- 
eled GW chirps with an interferometric detector network. 
Hence, we give more emphasis on the basic formalism and 
keep the colored noise case with basic minimal assump- 
tion: the noise from the different detectors is colored but 
with same covariance. Based on this ground work, the 
work is in progress to extend this to the colored noise 
case with different noise covariances. 

Let us therefore assume now that the noise compo- 
nents in each detector are independent and colored, 
with the same covariance matrix Rq. Recall that the 
covariance matrix of a random vector a is defined as 



The Levi-Civita symbol is defined as 

£ij... = +1 when i, j, . . . is an even permutation of 1, 2, 
= — 1 when i, j, . . . is an odd permutation of 1, 2, . 
= when any two labels are equal. 



E [(a - E[a])(a - E[a]) H ] where E[.] denotes the expec- 
tation. From the independence of the noise components, 
the overall covariance matrix of the network noise vector 
is then a block-diagonal matrix, where all the blocks are 
identical and equal to Ro: R = diag(Ro) = I<j ® Ro- 
In this case, the network LLR becomes 



A(x) 



where the notation 



np 



R- 



iR-i' 



(64) 



IR- 1 



denotes the norm induced 
by the inner product associated to the covariance matrix 
R _1 , i.e., Hajl^ = a ff R _1 a. 

Introducing the whitened version n = R 1,/2 II and 
x = R~ 1//2 x of II and x respectively, Eq. (64 1 can be 
rewritten as 



A(x) 



np 



■ 1 1 S3 

x . 



(65) 



which is similar to Eq. ( 38 1 where all the quantities are 



replaced by their whitened version. Thus, the maximiza- 
tion of A(x) with respect to the extrinsic parameters p 
will follow the same algebra as that derived in V How- 
ever, for the sake of completeness, we detail it below. 

Following Sec. [VJ maximizing A(x) with respect to the 
extrinsic parameters p leads to 



ri # i 



(66) 



where fi* is the pseudo-inverse of II. Expressing this 
pseudo-inverse by means of the SVD of II as II# = 
VjjS^U? and introducing Eq. ( 66 1 into Eq. (65l pro- 
vides the new statistic 



A(x) = ||Ufx|| 2 . 



(67) 



Now, from the definition of II and the specific structure 
of R, it is straightforward to see that 

ri = R~ 1/2 (D <g> E) = D ® (Rg 1/2 E) = D <g> E, (68) 
where we have introduced the whitened version e = 



R ( 



-1/2 



e of the chirp signal and the corresponding matrix 



E = [e e*]/2. 



For the white noise case, the statistic (67 1 can then be 



rewritten in terms of the SVD of the matrices D and E: 

' ~ 1 1 2 



A(x)= ||(U D (8)U^) ff x 



(69) 



The computation of Ug is similar to the computa- 
tion of U^. Furthermore, if we note that e T e ~ 0, 
and if we assume that e H e = N, it turns out that 
= 2E/VN = 2R~ 1/2 E,/y/N. Using the property 
of the Kronecker product in Eq. (A2|, we then obtain 



A(x) = ||vec(E ff XU*b)| 



(70) 



where the matrix X = [X1.X2 • • • x^] contains the data vec- 
tor from each detector whitened twice: x^- = R _1 / 2 Xj = 
R" x x,. 
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As this expression is similar to the white noise case, 

we can form two synthetic streams [yi,y2] = XU^, and 
use them to express the LLR statistic as 



A(x) = -(|e ff yi | 2 + |e^y 2 | 2 ). (71) 



In this expression, the only difference with the white 



noise LLR of (53 1 comes from the computation of the 



synthetic streams yi and yi which are obtained after 
double-whitening the data. 



VI. MAXIMIZATION OVER THE INTRINSIC 
PARAMETERS 

In the previous section, we maximized the network 
LLR over the extrinsic parameters of the signal model, 
assuming that the remaining parameters (the source lo- 
cation angles <f> and 9 and the phase function <^>(-)) were 
known. 

By definition, the intrinsic parameters modify the net- 
work LLR non-linearly For this reason, the maximiza- 
tion of A over these parameters is more difficult. It can- 
not be done analytically and must be performed numeri- 
cally, for instance with an exhaustive search of the max- 
imum by repeatedly computing A over the entire range 
of possibilities. 

While the exhaustive search can be employed for the 
source location angles, it is not applicable to the chirp 
phase function, which requires a specific method. For 
the single detector case, we had addressed this issue in 
4 with an original maximization scheme which is the 
cornerstone of the best chirplet chain (BCC) algorithm. 
Here, we use and adapt the principles of BCC to the 
multiple detector case. 



A. Chirp phase function 

Let us examine first the case of the detection of in- 
spiralling binary chirps. In this case, the chirp phase is 
a prescribed function of a small number of parameters 
i.e., the masses and spins of the binary stars. The maxi- 
mization over those is performed by constructing a grid of 
reference or template waveforms which are used to search 
the data. This grid samples the range of the physical pa- 
rameters. This sampling must be accurate (the template 
grid must be tight) to avoid missing any chirp. 

Tight grid of templates can be obtained in the non- 
parametric case (large number of parameters) i.e., when 
the chirp is not completely known. We have shown in [4 
how to construct a template grid which covers entirely 
the set of smooth chirps i.e., chirps whose frequency evo- 
lution has some regularity as described in Sec. |IIB 2| In 
the next section, we briefly describe this construction. 



1. Chirplet chains: tight template bank for smooth chirps 

We refer to the template forming this grid as chirplet 
chain (CC). These CCs are constructed on a simple geo- 
metrical idea: a broken line is a good approximation of a 
smooth curve. Since the frequency of a smooth chirp fol- 
lows a smooth frequency vs. time curve, we construct 
templates that are broken lines in the time-frequency 
(TF) plane. 

More precisely, CC are defined as follows. We start by 
sampling the TF plane with a regular grid consisting of 
N t time bins and Nf frequency bins. We build the tem- 
plate waveforms like a puzzle by assembling small chirp 
pieces which we refer to as chirplets. A chirplet is a signal 
with a frequency joining linearly two neighboring vertices 
of the grid. The result of this assembly is a chirplet chain 
i.e., a piecewise linear chirp. Since we are concerned with 
continuous frequency evolution with bounded variations, 
we only form continuous chains. 

We control the variations of the CC frequency. The fre- 
quency of a single chirplet does not increase or decrease 
more than N' r frequency bins over a time bin. Similarly, 
the difference of the frequency variations of two succes- 
sive chirplets in a chain does not increase or decrease 
more than TV" frequency bins. 

The CC grid is defined by four parameters namely Nt , 
Nf, N' r and N". Those are the available degrees of free- 
dom we can tune to make the CC grid tight. A tem- 
plate grid is tight if the network ambiguity A(s; <p') cx 
|e' yi(s)| 2 + |e' y2(s)| 2 which measures the similarity 
between an arbitrary chirp (of phase tp) and its closest 
template (of phase tp'), is large enough and relatively 
closed to the maximum (when tp' — tp). 

As stated in Eq. ( |55| . in presence of a noise free GW 
chirp, the synthetic streams y;(s) are rescaled and phase 
shifted copies of the initial GW chirp s. Therefore, we 
treat the network ambiguity as a sum of ambiguities from 
two virtual detectors where each terms in A(s; tp') is the 
ambiguity computed in the single detector case. An es- 
timate of the ambiguity has been obtained in [4] for this 
case. It can thus be directly reused to compute A(s; tp'). 

The bottom line is that the ratio of the ambiguity to 
its maximum for the network case remains unaltered as 
compared to the single detector case and thus same for 
the tight grid conditions. In conclusion, the rules (which 
we won't repeat here) established in [3] to set the search 
parameters can also be applied here. 



2. Search through CCs in the time-frequency plane: best 
network CC algorithm 

We have now to search through the CC grid to 
find the best matching template, i.e., which maximizes 
max all CCs v ' A(x; tp'). Counting the number of possible 
CCs to be searched over is a combinatorial problem. This 
count grows exponentially with the number of time bins 
Nt. In the situation of interest, it reaches prohibitively 
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large values. The family of CCs cannot be scanned ex- 
haustively and the template based search is intractable. 

In [4J , we propose an alternative scheme yielding a close 
approximation of the maximum for the single detector 
statistic. When applied to the network, the scheme de- 
mands to reformulate the network statistic in the TF 
plane. The TF plane offers a natural and geometrically 
simple representation of chirp signals which simplifies the 
statistic. It turns out that the resulting statistic falls in a 
class of objective functions where efficient combinatorial 
optimization algorithms can be used. We now explain 
this result in more detail. 

We use the TF representation given by the discrete 
Wigner-Ville (WV) distribution [27 defined for the time 
series x[n] with n = 0, . . . , JV — 1 as 

+fc„ 

w x (n,m)= x[\n+k/2\]x*[\n~k/2\]e- 2 " mk/{2N \ 

(72) 

with k n = min{2n,27V — 1 — 2n}, where |_-J gives the 
integer part. The arguments of w x are the time index n 
and the frequency index m which correspond, in physical 
units, to the time t n — t s n and the frequency v m — 
v 8 m/{2N) for < m < N and u m = u s (N - m)/(2N) 
for JV+l<m<2iV-l. 

The above WV distribution is a unitary representation. 
This means that the scalar products of two signals can 
be re-expressed as scalar products of their WV. Let x\[n] 
and X2[n] be two time series. The unitarity property of 
w x is expressed by the Moyal's formula as stated below 



Thus, substituting in Eq. (74 1, we obtain the following 



N-i 

n=0 



I N-12N-1 

-T7 ^2 ^2 w Xl (n,m)w X2 (n,Tn). 



2N 



n—0 rn—0 



(73) 

Applying this property to the network statistic in 
Eq. ( 53 ) , we get 



reformulation of the network statistic 



N-l 



(76) 



n=0 



The maximization of A(x) over the set of CC amounts 
to finding the TF path that maximizes the integral 
Eq. (76 1, which is equivalent to a longest path problem 



in the TF plane. This problem is structurally identical 
to the single detector case (the only change is the way we 
obtain the TF map) . We can therefore essentially re-use 
the scheme proposed earlier for this latter case. The lat- 
ter belongs to a class of combinatorial optimization prob- 
lems where efficient (polynomial time) algorithms exist. 
We use one such algorithm, namely, the dynamic pro- 
gramming. 

In conclusion, the combination of the two ingredients 
namely the synthetic streams and the phase maximization 
scheme used in BCC allows us to coherently search the 
unmodeled GW chirps in the data of GW detector net- 
work. We refer to this procedure as the best network CC 
(BNCC) algorithm. 



B. Source sky position 

As we are performing maximization successively, till 
now we assume that we know the sky position of the 
source. Knowing the sky position, we construct the 
synthetic streams with appropriate direction dependent 
weight factors, time-delay shifts and carry out the BNCC 
algorithm for chirp phase detection. In reality, the sky 
position is unknown. One needs to search through the 
entire sky by sampling the celestial sphere with a grid 
and repeating the above procedure for each point on this 
grid. 



^ N-12N-1 

^W = JpYl H w y (n,m)w e (n,m). 

n—0 m—0 



(74) 



where w y — w yi + w V2 combines the individual WVs of 
the two synthetic streams. 

In order to compute A(x), we need to have a model for 
w e . We know that the WV distribution of a linear chirp 
(whose frequency is a linear function of time) is essen- 
tially concentrated in the neighborhood of its instanta- 
neous frequency [27] . We assume that it also holds true 
for an arbitrary (non-linear) chirp. Applying this approx- 
imation to the WV w e of the template CC in Eq. (74 1, 
we get 



w e (n, to) « 2N S(m — m n ), 



(75) 



where m„ denotes the nearest integer of 2T v{t n ) and v 
is the instantaneous frequency of the CC. 



C. Time of arrival 

Since we process the data streams sequentially and 
block- wise, the maximization over to amounts to selecting 
that block where the statistic arrives at a local maximum 
(i.e., the maximum of the "detection peak"). The epoch 
of this block yields an estimate of to. The resolution of 
the estimate may be improved by increasing the overlap 
between two consecutive blocks. 



D. Estimation of computational cost 

We estimate the computational cost of the BNCC 
search by counting the floating-point operations (flops) 
required by its various subparts. The algorithm con- 
sists essentially in repeating the one-detector search for 
all sky location angles. Let Nq be the number of bins 
of the sky grid. The total cost is therefore Nq times 
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the cost of the one-detector search, which we give in jl] 
and summarize now. The computation of the WV of the 
two synthetic streams requires lONNf log 2 Nf flops and 
the BCC search applied to the combined WV requires 
[N + (2N'J + l)N t ]N c flops, where N c m {2N' r + l)N f is 
the total number of chirplets. Since this last part of the 
algorithm dominates, the over-all cost thus scales with 

C oc N n [N + + l)N t ](2N' r + l)N f . (77) 

These operations are applied to each blocks (of du- 
ration T) the data streams. The computational power 
needed to process the data in real time is thus given by 
C/(pT) where fi is the overlap between two successive 
blocks. 



VII. RESULTS WITH SIMULATED DATA AND 
DISCUSSION 

A. Proof of principle of a full blind search 

We present here a proof of principle for the proposed 
detection method. For this case study, we consider a 
network of three detectors placed and oriented like the 
existing Virgo and the two four-kilometer LIGO detec- 
tors. The coordinates and orientation of these detectors 
can be found in Table |TJ We assume a simplified model for 
detector noise which we generate independently for each 
detector, using a white Gaussian noise. Fig. [5] illustrates 
the possibility of a "full blind" search in this situation. 
This means that we perform the detection jointly with 
the estimation of the GW chirp frequency and the source 
sky location. 

1. Description of the test signal 

Because of computational limitation, we restrict this 
study to rather short chirps of N = 256 samples, i.e., 
a chirp duration T = 250 ms assuming a sampling rate 
of v s = 1024 Hz. The chirp frequency follows a random 
time evolution which however satisfies chirping rate con- 
straints. We make sure that the first and second deriva- 
tives of the chirp frequency are not larger than F' = 9.2 
kHz/s and F" = 1.57 MHz/s 2 respectively. The chosen 
test signal has about 50 cycles. This is a larger number 
than what is considered typically for burst GWs ( ~ 10). 

As a comparison with a well-known physical case, an 
inspiralling (equal mass) binary with total mass M = 
IIMq reaches the same maximum frequency variations at 
the last stable circular orbit. (Binary chirps with larger 
total mass also satisfy these chirping rate limits). 

We set the SNR to p = 20. The chirp is injected at 
the sky position <p = 2.8 rad and 9 — 0.4 rad where 
the contributions of the individual detectors are compa- 
rable, namely the individual SNR are 10.4, 10.15 and 
13.77 for Virgo, LIGO Hanford and LIGO Livingstone 
respectively. 



Ikelihcod tendErapfl injedcn poirtl +: 
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FIG. 5: Color online — Coherent detec- 

tion/estimation of a "random" chirp with network of 
three GW antennas. In the data of three GW antennas 
(the two LIGO and Virgo), we inject (a) a "random" GW 
chirp emitted from a source at the position marked with 
a "+" at tf> = 2.8 rad and 6 = 0.4 rad. We perform a full 
sky search using the best network chirplet chain algorithm. 
It produces a likelihood landscape (b) where we select the 
maximum. This is the detection point and it is indicated 
with "x". In (c) we show the combined WV distribution 
of the synthetic streams at the detection point. In (d), we 
compare the exact frequency of the chirp (solid/blue) with 
the estimation (dashed/red) obtained at the detection point. 



2. Search parameters 



We search through the set of CCs defined over a TF 
grid with N t = 128 time intervals and N = Nf = 256 
frequency bins (using f s = 2048 Hz). We set the regular- 
ity parameters to Nf. = 9 and N" = 3, consistent to the 
above chirping rate limits. 

We select an ad-hoc sky grid by dividing regularly the 
full range of the source localization angles 9 and <f> into 
128 bins. The resulting grid has therefore a total of 
Nq = 16384 bins. This is probably much finer than is 
required to perform the detection without missing can- 
didate. However, this oversampling leads to precise like- 
lihood sky maps which helps to diagnose the method. 
With this parameter choice, the estimated computational 
power required to analyzed the data in real time is of 2.8 
TFlops, assuming an overlap of fi = 50% between suc- 
cessive data blocks. Because of the crude choice for the 
sky grid, this requirement is probably over-estimated. 

The result of the search is displayed in Fig. [5] where 
we see that the injection is recovered both in sky po- 
sition and frequency evolution. The source position is 
estimated at cj> = 2.8 rad and 9 = 0.39 rad. 
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B. Regularized variants 



likelitood fcandscspe, injedbn poinl +: siimaled Iraq, al dsledkm pdm x: 

4=0 71 rad, 6=233 rad, corxJ=192 61 t=0S7 rad ft=2.63 rad 

1000 F 



As shown in Sec. V B 1| the SNR carried by the syn- 
thetic stream is proportional to the corresponding SV. 
When the GW polarization plane is degenerate (i.e., 
when (72 is small), the second synthetic stream contains 
almost only noise. We thus don't lose information if we 
suppress its contribution from the statistic. This is the 
basic idea of Klimenko et al. in [9]. 

We have seen that the estimation of the extrinsic pa- 
rameters is an ill-posed least-square problem in those 
cases. Suppressing the contribution of the second syn- 
thetic stream amounts to regularizing this problem [5]. 
In practice, this regularization can be done in various 
ways, corresponding to well-identified schemes. 

A first possibility is to suppress the contribution of 
the second synthetic stream when the conditioning num- 
ber of II is too large (i.e., exceeds a given threshold). 
This scheme is referred to as truncated SVD [23]. A sec- 
ond possibility is to balance (divide) the contribution of 
the second synthetic stream by the conditioning number. 
This is referred to as the Tykhonov approach [23J and it 
was proposed for regularizing burst searches in [TT] , 

In Fig. [6) we compare the likelihood landscape and fre- 
quency estimate obtained with the standard statistic and 
its regularized version using the Tykhonov approach. Vi- 
sually, the regularization improves the contrast and con- 
centration of the likelihood landscape around the injec- 
tion point. This can be assessed more quantitatively with 
the contrast defined as the ratio of the likelihood land- 
scape extremes. This contrast is improved by about 10% 
for the regularized statistic as compared to the standard 
version. It is also interesting to compare the "width" of 
the detection peak obtained with the two statistic. To 
do this, we measure the solid angle of the the sky region 
where the statistic is larger than 90% of the maximum. 
This angle is reduced by a factor of ~ 6 when computed 
with the regularized statistic. There is however no major 
improvement of the frequency estimate. More generally, 
it is unclear whether the regularized statistic performs 
better than the standard one. 



VIII. CONCLUDING REMARKS 

The coherent detection of unmodclcd chirps with a net- 
work of GW have features and issues in common with 
the burst one. In particular, the same geometrical ob- 
jects play a key role. While the noise spans the whole 
d dimensional data space, GW signals (chirps or burst) 
only belong to a two-dimensional (one dimension per GW 
polarizations) subspace, the GW polarization plane. De- 
tecting GWs amounts to checking whether the data has 
significant components in this plane or not. To do so, 
we compute the projections of the data onto a basis of 
the GW polarization plane. In practice, this defines two 
instantaneous linear mixtures of the individual detector 
data which we refer to as synthetic streams. Those may 




FIG. 6: Color online — Coherent detec- 

tion/estimation of a "random" chirp with network 
of three GW antennas. Standard and regularized 
statistic. We compare the likelihood landscape (left) and 
frequency estimation (right) obtained using the standard 
(top) and Tykhonov-regularized (bottom) versions of the 
network statistic. The test signal is a "random" GW chirp 
injected at the sky location marked with a "+". This location 
has been chosen because of the associated large value of the 
conditioning number, namely cond (II) w 192. 



be considered as the output data of "virtual" detectors. 
This combination is such that the GW contributions from 
each real detector add constructively. The GW signature 
thus has a larger amplitude in the synthetic streams while 
the noise variance is kept at the same level. 

The coherent detection amounts to looking for an ex- 
cess in the signal energy in one or both synthetic streams 
(depending on the GW polarization model). This pro- 
vides a generic and simple procedure to produce a coher- 
ent detection pipeline from a one detector pipeline. In 
the one detector case, the BCC search performs a path 
search in a time-frequency distribution of the data. In 
the multiple detector case, the BNCC search now use 
the joint time-frequency map obtained by summing the 
time-frequency energy distributions of the two synthetic 
streams. The approach does not restrict to chirp detec- 
tion and it can be applied to burst searches [14]. 

We demonstrated in a simplified situation that the full- 
sky blind detection of an unmodcled chirp is feasible. 
This means that the detection is performed jointly to 
estimate the source location and the frequency evolution. 
The application of this method to the real data however 
requires several improvements. First, the method has to 
be adapted to the case where the detectors have different 
sensitivities. In this respect, we already obtained first 
results [28] . 

We also have to refine the choice of the grid which sam- 
ples the celestial sphere. In the present work, we select 
an ad-hoc grid. Clearly, the sky resolution and the bin 
shape depends on the geometry of the sphere and the 
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location and orientation of the detectors in the consid- 
ered network. A better grid choice (not too coarse to 
avoid SNR loses, or and not too fine to avoid using use- 
less computing resources) should incorporate this infor- 
mation keeping the search performance (detection prob- 
ability and sky resolution) constant. In this respect, we 
may consider to other parameterizations of the sky lo- 
cation which makes the definition of the sky grid easier, 
for instance by choosing the time-delays as investigated 
in [7]. We may also explore hierarchical schemes for the 
reduction of the computational cost. 

The GW polarization plane depends on the detector 
antenna patterns functions. With the presently available 
networks, there are significantly large sky regions where 
the antenna patterns are almost aligned. In this case, the 
network observes essentially one polarization and is al- 
most insensitive to the other: the GW polarization plane 
reduces to a one-dimensional space. The information car- 
ried by the missing polarization lacks and this makes the 
estimation of certain parameters ill-posed and hence very 
sensitive to noise. We can evaluate that the variance of 
the estimate scales with the condition number of the an- 
tenna pattern matrix. When this number (which quanti- 
fies in some sense the mutual alignment of the detectors) 
is large, the estimation is ill-posed and we expect poor 
results. 

This is an important issue for burst detection since it 
affects significantly the shape of the estimated waveform 
(and, particularly the regularity of its time evolution). 
This has motivated the development of regularization 
schemes which penalize the estimation of non-physical 
(i.e., irregular) waveforms. We have shown that this is 
however less of a problem for chirps because of their more 
constrained model. Ill-conditioning only affects global 
scaling factors in the chirp model. Unlike bursts, no ad- 
ditional prior is available for regularizing the estimation 
of these scaling factors. 

The data space can be decomposed as the direct sum 
of the GW polarization plane and its complementary. 
While GWs have zero components in the latter null space, 
it is unlikely that instrumental noise (including its non- 
Gaussian and non-stationary part) will. This motivates 
the use of null streams (i.e., the projection of the data 
along a basis of the complementary space) to verify that 
a trigger is indeed a GW candidate and not an instru- 
mental artifact. Since null streams are inexpensive to 
compute, we consider to use them to make preemptive 
cuts in order to avoid the analysis of bad data. 
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APPENDIX A: KRONECKER PRODUCT: 
DEFINITION AND PROPERTIES 



The Kronecker product ® transforms two matrices A € 
: mxn and B <E C pxq into the following matrix of C mpxnq 



A (g) B = 



anB 



a m iB 



a ln B 



(Al) 



The Kronecker product is a linear transform and can 
be considered as a special case of the tensor product. 
We define the operator vec() to be the stack operator 
which transforms the matrix into a vector by stacking 
its columns, i.e. x = vec(X). In the text, we use the 
following property: 

(A ® B)vec(X) = vec(BXA T ) . (A2) 

The proof of this property is straightforward. 



APPENDIX B: INTERFEROMETRIC 
DETECTOR RESPONSE IN TERMS OF 
GEL'FAND FUNCTIONS 

The GW response of a detector to an incoming GW 
can be obtained by computing the interaction of the wave 
tensor W with the detector tensor D as follows 12 : 



3 



(Bl) 



The wave tensor is related to the incoming GW tensor 
in the TT gauge by — 2Wy . Both detector and wave 
tensors are rank 2 STF tensors. Any STF tensor can be 
expanded in the basis of spin- weighted spherical harmon- 
ics of rank 2 and the rank-2 Gel'fand functions provide 
the corresponding coefficients. Further, they are repre- 
sentation of rotation group SO (3) and provide compact 
representation for the detector response of any arbitrarily 
oriented and located detector on Earth which we present 
in this appendix |29j . 

a. Wave tensor — The incoming GW tensor in TT 
gauge is given by 

h>ij — {exi^xj ^yi^yj)ti-\- ~t~ 2(e a; zGyj)/Z'x , (-^^) 

where e x and e y are unit vectors along the x w and y w 
axes in the wave frame; h + and h x are the two GW 
polarizations. Let to = (e x + ie y )/^/2 be a complex 



12 Unless otherwise mentioned, the notations and symbols used in 
all appendices are confined to those appendices only. 
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vector in the wave frame. Then, the wave tensor can be 
written down in terms of m as 



(B3) 
■ih x which 



where we used the complex quantity h = h-\ 
combines both GW polarizations. 

The term mirrij is a STF tensor of rank 2. We choose to 
work in the detector frame for convenience. Expanding 
niirrij in terms of spin- weighted spherical harmonics of 
rank 2, namely y 2n and the rank-2 Gel'fand functions 
T mn {(f>' e ,0' e ,ip' e ), we get 



2 

:r> ^ 

n= — 2 



(B4) 



The angles ((t>' c . 6{ , ip' e ) are the Euler angles of the rotation 
operator which transforms the detector coordinates into 
the wave coordinates. 



Substituting Eq. (B4| into Eq. (B3l, we express the 



wave tensor in terms of the Gel'fand functions as 

2 



W 



nwtT-in] ■ 

n=-2 

The detector tensor is 



b. Detector tensor 

Dij = n u nij - n 2l n 2j , 



(B5) 



(B6) 



where ni and n 2 are the unit vectors along the first and 
second arms of the interferometer. Recall that we choose 
the Xd-axis of the detector frame along the bisector of the 
two arms. The y^-axis is chosen such that (xd, yd, Zd) is a 
right handed coordinate system with Zd pointing towards 
the local zenith. With this choice, we have 



D 22 = 0, D 12 = D 21 = -l. (B7) 

From Eqs. (B6l and ( |B7[ ), the detector response is 
s = $t[f*h], 



(B8) 



where the complex antenna pattern function is given by 

/ = i[T 2 _ 2 (0' e ,CVe) " TMXrfe)] ■ (B9) 



From the expansion of above Eq. ( B8 1 in terms of GW 
polarizations, it is consistent to define f = f+ + ifx, 
which yields Eq. ( 13 1 . 



c. Extended complex antenna pattern for sources or- 
biting in a fixed plane — As discussed in Sec. |IIID 1[ 
the extended antenna pattern functions incorporate the 



inclination angle e and is given in Eq. ( 17 1 as 



- l + cos 2 e . 

/ = o /+ + 1 cos e /x 



(BIO) 



where /+ jX depend on the relative orientation of wave 
frame with respect to the detector frame. The detector 
to wave frame coordinate transformation can be split into 



two: detector to Earth's frame and Earth's frame to wave 
frame by the following rotation transformations as given 
in Eqs. Q and Q 

e' e , i>' e ) - o(</> e , e e , yjo-^a, 7 ) • (bii) 

The above successive rotation transformation can be 
translated into the addition theorem of Gel'fand func- 
tions [30 as given below 



Z=-2 

(B12) 

We used the fact that the inverse rotation operator is 
associated to a complex conjugation. 



Substituting in Eq. ( B9 ) , we rewrite the antenna pat- 



tern functions in terms of the Gel'fand functions as 

2 

/= iT 2s(<l)e,de,tpe)[T 2s (a,P,j) - T_ 2s (a, /3, 7)]* . 

s=-2 

_ (B13) 

Substituting in the extended beam pattern function 
given in Eq. (BIO ) and combining the dependencies upon 



tj) = ijj e and e, we get 

/ = T 22 (i>, e, Q)d + T 2 - 2 (iP, e, 0)d* , (B14) 

where 

2 

d = - iT 2n (<f) e ,6 e ,0)[T 2n (a,j3,j) - T_ 2n (a,^,7)]* . 
n=— 2 

(B15) 

There are various ways of expressing the antenna pat- 
tern functions. The main advantages of this one is that 
it is particularly compact and that the angles "0 and e get 
factorized from the rest of the parameters. This helps in 
the maximization of the network LLR over the extrinsic 
parameters. 



APPENDIX C: SVD OF A TWO-COLUMN 
COMPLEX MATRIX 

In this appendix, we obtain the SVD of a complex ma- 
trix of the type A = [a, a*] where a E C Nxl . The SVD 
decomposes A into the product A = U a^a^a where 
U,4 and V a are two orthogonal matrices and Ha is a pos- 
itive definite diagonal matrix. To obtain it analytically, 
we first get the eigen-decomposition of 



A H A = 





a ff a*" 




a 


b* 


T 

a a 


T * 

a a 




b 


a 



(CI) 



We distinguish two cases depending on the value of b. 
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1. For |fc| > 

The eigenvalues a± and the eigenvectors v± of A H A 
are given below; 



4 = o±|6|, 
1 

71 



exp(— ivj) 
±1 



(C2) 
(C3) 



where vo = arg& = arg[a a]. The number of non-zero 
eigenvalues of A H A gives the rank of A. The eigenval- 
ues arranged in the descending order form the diagonal 
values of Y? A as shown below. The eigenvectors v± are 
the right-handed singular vectors of A and form V^4 as 
given below: 



S A = 

v A = 



(7 _|_ 

a- 
v + v 



(C4) 
(C5) 



We then form the matrix = AV^Sj = [u + , u_] 
containing the right-handed singular vectors, with 



u± 



exp(— izu)a ± a* 
V2<J± 



(C6) 



The above expressions are valid only when the two 
SV are non zero. It is possible that the smallest SV cr_ 
vanishes. In this case, the SVD collapses to XJa = u + , 
E A = <j + , and V A = v + . 



For 161 = 



In this case, the matrix A H A = di.2 is diagonal. We 
thus have £ = s/a^ , "Va = 1 2 and — A/ y/a. 
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